﻿///////////////////////////////////////////////////////////////////////////////
//
//  This file is part of MathLib.NET.
//
//  This library is free software; you can redistribute it and/or
//  modify it under the terms of the GNU Lesser General Public
//  License as published by the Free Software Foundation; either
//  version 2.1 of the License, or (at your option) any later version.
//  
//  This library is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//  Lesser General Public License for more details.
//  
//  You should have received a copy of the GNU Lesser General Public
//  License along with this library;  If not, see 
//  <http://www.gnu.org/licenses/>.
//
///////////////////////////////////////////////////////////////////////////////

using System.Collections.Generic;
using ILNumerics;
using ILNumerics.BuiltInFunctions;

namespace MathLib
{
    /// <summary>
    /// Helper class for accessing ILNumeric arrays.
    /// </summary>
    public class ArrayHelper
    {
        /// <summary>
        /// Given an array and leading dimension, this method returns the number of transforms required,
        /// the number of elements between the start of adjacent transforms in the internal array, the
        /// number of points in each transform and number of elements between adjacent transform points.
        /// </summary>
        /// <typeparam name="T">Array cell type</typeparam>
        /// <param name="array">Array to calculate parameters for</param>
        /// <param name="dim">Leading dimension along which to perform transformations</param>
        /// <param name="numTransforms">Total number of transforms to perform</param>
        /// <param name="majorTransformStride">Number of elements from one transform to the next for dimensions more significant than dim</param>
        /// <param name="numPoints">Number of points per transform</param>
        /// <param name="pointStride">Number of elements from one point to the next</param>
        public static void CalcArrayAccessParams<T>(ILArray<T> array, ref int dim, out int numTransforms, out int majorTransformStride, out int numPoints, out int pointStride)
        {
            // if dim = -1 use first non-singleton dimension
            if (dim == -1)
                dim = array.Dimensions.FirstNonSingleton();

            // If dim = -1 (scalar array) use dimension 1
            if (dim == -1)
                dim = 1;

            // Number of points = length of supplied dimension
            numPoints = array.Dimensions[dim];

            // Number of transforms = total number of elements / number of points
            numTransforms = array.Dimensions.NumberOfElements / numPoints;

            // Point stride = product of dimensions up to but not including dim
            pointStride = 1;
            for (int i = 0; i < dim; i++)
                pointStride *= array.Dimensions[i];

            // Transform stride - product of dimensions including dim
            majorTransformStride = pointStride * array.Dimensions[dim];
        }

        /// <summary>
        /// Performs cross multiplication of the two supplied real arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Multiplication of two arrays</returns>
        public static ILArray<double> multiply(ILArray<double> a1, ILArray<double> a2)
        {
            ILArray<double> output;

            if (a1.IsScalar)
            {
                // Scalar product
                output = new ILArray<double>(a2.Dimensions.ToIntArray());
                for (int i = 0; i < a2.InternalArray4Experts.Length; i++)
                {
                    output.InternalArray4Experts[i] = a1.InternalArray4Experts[0] * a2.InternalArray4Experts[i];
                }
            }
            else if (a2.IsScalar)
            {
                // Scalar product
                output = new ILArray<double>(a1.Dimensions.ToIntArray());
                for (int i = 0; i < a1.InternalArray4Experts.Length; i++)
                {
                    output.InternalArray4Experts[i] = a2.InternalArray4Experts[0] * a1.InternalArray4Experts[i];
                }
            }
            else if (a1.Dimensions.NumberOfDimensions == 2 && a2.Dimensions.NumberOfDimensions == 2)
            {
                // Matrix multiplication
                int a1c = a1.Dimensions[1];
                int a1r = a1.Dimensions[0];
                int a2c = a2.Dimensions[1];
                int a2r = a2.Dimensions[0];
                if (a1c == a2r)
                {
                    output = new ILArray<double>(new int[] { a2c, a1r });
                    for (int r = 0; r < a1r; r++)
                        for (int c = 0; c < a2c; c++)
                            for (int i = 0; i < a1c; i++)
                                output.InternalArray4Experts[c * a1r + r] += a1.InternalArray4Experts[i * a1r + r] * a2.InternalArray4Experts[c * a2r + i];
                }
                else
                    throw new MathException("multiply(): matrix multiplication requires complementary matrices (i, k), (k, j).");
            }
            else
                throw new MathException("multiply(): requires scalar, vector or matrix parameters.");

            return output;
        }

        /// <summary>
        /// Performs cross multiplication of the two supplied real arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Multiplication of two arrays</returns>
        public static ILArray<float> multiply(ILArray<float> a1, ILArray<float> a2)
        {
            ILArray<float> output;

            if (a1.IsScalar)
            {
                // Scalar product
                output = new ILArray<float>(a2.Dimensions.ToIntArray());
                for (int i = 0; i < a2.InternalArray4Experts.Length; i++)
                {
                    output.InternalArray4Experts[i] = a1.InternalArray4Experts[0] * a2.InternalArray4Experts[i];
                }
            }
            else if (a2.IsScalar)
            {
                // Scalar product
                output = new ILArray<float>(a1.Dimensions.ToIntArray());
                for (int i = 0; i < a1.InternalArray4Experts.Length; i++)
                {
                    output.InternalArray4Experts[i] = a2.InternalArray4Experts[0] * a1.InternalArray4Experts[i];
                }
            }
            else if (a1.Dimensions.NumberOfDimensions == 2 && a2.Dimensions.NumberOfDimensions == 2)
            {
                // Matrix multiplication
                int a1c = a1.Dimensions[1];
                int a1r = a1.Dimensions[0];
                int a2c = a2.Dimensions[1];
                int a2r = a2.Dimensions[0];
                if (a1c == a2r)
                {
                    output = new ILArray<float>(new int[] { a2c, a1r });
                    for (int r = 0; r < a1r; r++)
                        for (int c = 0; c < a2c; c++)
                            for (int i = 0; i < a1c; i++)
                                output.InternalArray4Experts[c * a1r + r] += a1.InternalArray4Experts[i * a1r + r] * a2.InternalArray4Experts[c * a2r + i];
                }
                else
                    throw new MathException("multiply(): matrix multiplication requires complementary matrices (i, k), (k, j).");
            }
            else
                throw new MathException("multiply(): requires scalar, vector or matrix parameters.");

            return output;
        }

        /// <summary>
        /// Performs multiplication of the two supplied complex arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Multiplication of two arrays</returns>
        public static ILArray<complex> multiply(ILArray<complex> a1, ILArray<complex> a2)
        {
            ILArray<complex> output;

            if (a1.IsScalar)
            {
                // Scalar product
                output = new ILArray<complex>(a2.Dimensions.ToIntArray());
                for (int i = 0; i < a2.InternalArray4Experts.Length; i++)
                {
                    output.InternalArray4Experts[i] = a1.InternalArray4Experts[0] * a2.InternalArray4Experts[i];
                }
            }
            else if (a2.IsScalar)
            {
                // Scalar product
                output = new ILArray<complex>(a1.Dimensions.ToIntArray());
                for (int i = 0; i < a1.InternalArray4Experts.Length; i++)
                {
                    output.InternalArray4Experts[i] = a2.InternalArray4Experts[0] * a1.InternalArray4Experts[i];
                }
            }
            else
            {
                // Matrix multiplication
                int a1c = a1.Dimensions[1];
                int a1r = a1.Dimensions[0];
                int a2c = a2.Dimensions[1];
                int a2r = a2.Dimensions[0];
                if (a1c == a2r)
                {
                    output = new ILArray<complex>(new int[] { a2c, a1r });
                    for (int r = 0; r < a1r; r++)
                        for (int c = 0; c < a2c; c++)
                            for (int i = 0; i < a1c; i++)
                                output.InternalArray4Experts[c * a1r + r] += a1.InternalArray4Experts[i * a1r + r] * a2.InternalArray4Experts[c * a2r + i];
                }
                else
                    throw new MathException("multiply(): matrix multiplication requires complementary matrices (i, k), (k, j).");
            }

            return output;
        }

        /// <summary>
        /// Performs multiplication of the two supplied complex arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Multiplication of two arrays</returns>
        public static ILArray<fcomplex> multiply(ILArray<fcomplex> a1, ILArray<fcomplex> a2)
        {
            ILArray<fcomplex> output;

            if (a1.IsScalar)
            {
                // Scalar product
                output = new ILArray<fcomplex>(a2.Dimensions.ToIntArray());
                for (int i = 0; i < a2.InternalArray4Experts.Length; i++)
                {
                    output.InternalArray4Experts[i] = a1.InternalArray4Experts[0] * a2.InternalArray4Experts[i];
                }
            }
            else if (a2.IsScalar)
            {
                // Scalar product
                output = new ILArray<fcomplex>(a1.Dimensions.ToIntArray());
                for (int i = 0; i < a1.InternalArray4Experts.Length; i++)
                {
                    output.InternalArray4Experts[i] = a2.InternalArray4Experts[0] * a1.InternalArray4Experts[i];
                }
            }
            else
            {
                // Matrix multiplication
                int a1c = a1.Dimensions[1];
                int a1r = a1.Dimensions[0];
                int a2c = a2.Dimensions[1];
                int a2r = a2.Dimensions[0];
                if (a1c == a2r)
                {
                    output = new ILArray<fcomplex>(new int[] { a2c, a1r });
                    for (int r = 0; r < a1r; r++)
                        for (int c = 0; c < a2c; c++)
                            for (int i = 0; i < a1c; i++)
                                output.InternalArray4Experts[c * a1r + r] += a1.InternalArray4Experts[i * a1r + r] * a2.InternalArray4Experts[c * a2r + i];
                }
                else
                    throw new MathException("multiply(): matrix multiplication requires complementary matrices (i, k), (k, j).");
            }

            return output;
        }

        /// <summary>
        /// Performs multiplication of the two supplied complex arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Multiplication of two arrays</returns>
        public static ILArray<long> multiply(ILArray<long> a1, ILArray<long> a2)
        {
            ILArray<long> output;

            if (a1.IsScalar)
            {
                // Scalar product
                output = new ILArray<long>(a2.Dimensions.ToIntArray());
                for (int i = 0; i < a2.InternalArray4Experts.Length; i++)
                {
                    output.InternalArray4Experts[i] = a1.InternalArray4Experts[0] * a2.InternalArray4Experts[i];
                }
            }
            else if (a2.IsScalar)
            {
                // Scalar product
                output = new ILArray<long>(a1.Dimensions.ToIntArray());
                for (int i = 0; i < a1.InternalArray4Experts.Length; i++)
                {
                    output.InternalArray4Experts[i] = a2.InternalArray4Experts[0] * a1.InternalArray4Experts[i];
                }
            }
            else
            {
                // Matrix multiplication
                int a1c = a1.Dimensions[1];
                int a1r = a1.Dimensions[0];
                int a2c = a2.Dimensions[1];
                int a2r = a2.Dimensions[0];
                if (a1c == a2r)
                {
                    output = new ILArray<long>(new int[] { a2c, a1r });
                    for (int r = 0; r < a1r; r++)
                        for (int c = 0; c < a2c; c++)
                            for (int i = 0; i < a1c; i++)
                                output.InternalArray4Experts[c * a1r + r] += a1.InternalArray4Experts[i * a1r + r] * a2.InternalArray4Experts[c * a2r + i];
                }
                else
                    throw new MathException("multiply(): matrix multiplication requires complementary matrices (i, k), (k, j).");
            }

            return output;
        }

        /// <summary>
        /// Performs multiplication of the two supplied complex arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Multiplication of two arrays</returns>
        public static ILArray<ulong> multiply(ILArray<ulong> a1, ILArray<ulong> a2)
        {
            ILArray<ulong> output;

            if (a1.IsScalar)
            {
                // Scalar product
                output = new ILArray<ulong>(a2.Dimensions.ToIntArray());
                for (int i = 0; i < a2.InternalArray4Experts.Length; i++)
                {
                    output.InternalArray4Experts[i] = a1.InternalArray4Experts[0] * a2.InternalArray4Experts[i];
                }
            }
            else if (a2.IsScalar)
            {
                // Scalar product
                output = new ILArray<ulong>(a1.Dimensions.ToIntArray());
                for (int i = 0; i < a1.InternalArray4Experts.Length; i++)
                {
                    output.InternalArray4Experts[i] = a2.InternalArray4Experts[0] * a1.InternalArray4Experts[i];
                }
            }
            else
            {
                // Matrix multiplication
                int a1c = a1.Dimensions[1];
                int a1r = a1.Dimensions[0];
                int a2c = a2.Dimensions[1];
                int a2r = a2.Dimensions[0];
                if (a1c == a2r)
                {
                    output = new ILArray<ulong>(new int[] { a2c, a1r });
                    for (int r = 0; r < a1r; r++)
                        for (int c = 0; c < a2c; c++)
                            for (int i = 0; i < a1c; i++)
                                output.InternalArray4Experts[c * a1r + r] += a1.InternalArray4Experts[i * a1r + r] * a2.InternalArray4Experts[c * a2r + i];
                }
                else
                    throw new MathException("multiply(): matrix multiplication requires complementary matrices (i, k), (k, j).");
            }

            return output;
        }

        /// <summary>
        /// Performs multiplication of the two supplied complex arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Multiplication of two arrays</returns>
        public static ILArray<int> multiply(ILArray<int> a1, ILArray<int> a2)
        {
            ILArray<int> output;

            if (a1.IsScalar)
            {
                // Scalar product
                output = new ILArray<int>(a2.Dimensions.ToIntArray());
                for (int i = 0; i < a2.InternalArray4Experts.Length; i++)
                {
                    output.InternalArray4Experts[i] = a1.InternalArray4Experts[0] * a2.InternalArray4Experts[i];
                }
            }
            else if (a2.IsScalar)
            {
                // Scalar product
                output = new ILArray<int>(a1.Dimensions.ToIntArray());
                for (int i = 0; i < a1.InternalArray4Experts.Length; i++)
                {
                    output.InternalArray4Experts[i] = a2.InternalArray4Experts[0] * a1.InternalArray4Experts[i];
                }
            }
            else
            {
                // Matrix multiplication
                int a1c = a1.Dimensions[1];
                int a1r = a1.Dimensions[0];
                int a2c = a2.Dimensions[1];
                int a2r = a2.Dimensions[0];
                if (a1c == a2r)
                {
                    output = new ILArray<int>(new int[] { a2c, a1r });
                    for (int r = 0; r < a1r; r++)
                        for (int c = 0; c < a2c; c++)
                            for (int i = 0; i < a1c; i++)
                                output.InternalArray4Experts[c * a1r + r] += a1.InternalArray4Experts[i * a1r + r] * a2.InternalArray4Experts[c * a2r + i];
                }
                else
                    throw new MathException("multiply(): matrix multiplication requires complementary matrices (i, k), (k, j).");
            }

            return output;
        }

        /// <summary>
        /// Performs multiplication of the two supplied complex arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Multiplication of two arrays</returns>
        public static ILArray<uint> multiply(ILArray<uint> a1, ILArray<uint> a2)
        {
            ILArray<uint> output;

            if (a1.IsScalar)
            {
                // Scalar product
                output = new ILArray<uint>(a2.Dimensions.ToIntArray());
                for (int i = 0; i < a2.InternalArray4Experts.Length; i++)
                {
                    output.InternalArray4Experts[i] = a1.InternalArray4Experts[0] * a2.InternalArray4Experts[i];
                }
            }
            else if (a2.IsScalar)
            {
                // Scalar product
                output = new ILArray<uint>(a1.Dimensions.ToIntArray());
                for (int i = 0; i < a1.InternalArray4Experts.Length; i++)
                {
                    output.InternalArray4Experts[i] = a2.InternalArray4Experts[0] * a1.InternalArray4Experts[i];
                }
            }
            else
            {
                // Matrix multiplication
                int a1c = a1.Dimensions[1];
                int a1r = a1.Dimensions[0];
                int a2c = a2.Dimensions[1];
                int a2r = a2.Dimensions[0];
                if (a1c == a2r)
                {
                    output = new ILArray<uint>(new int[] { a2c, a1r });
                    for (int r = 0; r < a1r; r++)
                        for (int c = 0; c < a2c; c++)
                            for (int i = 0; i < a1c; i++)
                                output.InternalArray4Experts[c * a1r + r] += a1.InternalArray4Experts[i * a1r + r] * a2.InternalArray4Experts[c * a2r + i];
                }
                else
                    throw new MathException("multiply(): matrix multiplication requires complementary matrices (i, k), (k, j).");
            }

            return output;
        }

        /// <summary>
        /// Performs multiplication of the two supplied complex arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Multiplication of two arrays</returns>
        public static ILArray<short> multiply(ILArray<short> a1, ILArray<short> a2)
        {
            ILArray<short> output;

            if (a1.IsScalar)
            {
                // Scalar product
                output = new ILArray<short>(a2.Dimensions.ToIntArray());
                for (int i = 0; i < a2.InternalArray4Experts.Length; i++)
                {
                    output.InternalArray4Experts[i] = (short)(a1.InternalArray4Experts[0] * a2.InternalArray4Experts[i]);
                }
            }
            else if (a2.IsScalar)
            {
                // Scalar product
                output = new ILArray<short>(a1.Dimensions.ToIntArray());
                for (int i = 0; i < a1.InternalArray4Experts.Length; i++)
                {
                    output.InternalArray4Experts[i] = (short)(a2.InternalArray4Experts[0] * a1.InternalArray4Experts[i]);
                }
            }
            else
            {
                // Matrix multiplication
                int a1c = a1.Dimensions[1];
                int a1r = a1.Dimensions[0];
                int a2c = a2.Dimensions[1];
                int a2r = a2.Dimensions[0];
                if (a1c == a2r)
                {
                    output = new ILArray<short>(new int[] { a2c, a1r });
                    for (int r = 0; r < a1r; r++)
                        for (int c = 0; c < a2c; c++)
                            for (int i = 0; i < a1c; i++)
                                output.InternalArray4Experts[c * a1r + r] += (short)(a1.InternalArray4Experts[i * a1r + r] * a2.InternalArray4Experts[c * a2r + i]);
                }
                else
                    throw new MathException("multiply(): matrix multiplication requires complementary matrices (i, k), (k, j).");
            }

            return output;
        }

        /// <summary>
        /// Performs multiplication of the two supplied complex arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Multiplication of two arrays</returns>
        public static ILArray<ushort> multiply(ILArray<ushort> a1, ILArray<ushort> a2)
        {
            ILArray<ushort> output;

            if (a1.IsScalar)
            {
                // Scalar product
                output = new ILArray<ushort>(a2.Dimensions.ToIntArray());
                for (int i = 0; i < a2.InternalArray4Experts.Length; i++)
                {
                    output.InternalArray4Experts[i] = (ushort)(a1.InternalArray4Experts[0] * a2.InternalArray4Experts[i]);
                }
            }
            else if (a2.IsScalar)
            {
                // Scalar product
                output = new ILArray<ushort>(a1.Dimensions.ToIntArray());
                for (int i = 0; i < a1.InternalArray4Experts.Length; i++)
                {
                    output.InternalArray4Experts[i] = (ushort)(a2.InternalArray4Experts[0] * a1.InternalArray4Experts[i]);
                }
            }
            else
            {
                // Matrix multiplication
                int a1c = a1.Dimensions[1];
                int a1r = a1.Dimensions[0];
                int a2c = a2.Dimensions[1];
                int a2r = a2.Dimensions[0];
                if (a1c == a2r)
                {
                    output = new ILArray<ushort>(new int[] { a2c, a1r });
                    for (int r = 0; r < a1r; r++)
                        for (int c = 0; c < a2c; c++)
                            for (int i = 0; i < a1c; i++)
                                output.InternalArray4Experts[c * a1r + r] += (ushort)(a1.InternalArray4Experts[i * a1r + r] * a2.InternalArray4Experts[c * a2r + i]);
                }
                else
                    throw new MathException("multiply(): matrix multiplication requires complementary matrices (i, k), (k, j).");
            }

            return output;
        }

        /// <summary>
        /// Performs multiplication of the two supplied complex arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Multiplication of two arrays</returns>
        public static ILArray<double> multiply(ILArray<char> a1, ILArray<char> a2)
        {
            return multiply(ILMath.todouble(a1), ILMath.todouble(a2));
        }

        /// <summary>
        /// Performs multiplication of the two supplied complex arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Multiplication of two arrays</returns>
        public static ILArray<double> multiply(ILArray<byte> a1, ILArray<byte> a2)
        {
            return multiply(ILMath.todouble(a1), ILMath.todouble(a2));
        }

        /// <summary>
        /// Performs left division of the two supplied complex arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Left-division of two arrays</returns>
        public static ILArray<complex> ldivide(ILArray<complex> a1, ILArray<complex> a2)
        {
            // Convert matrices to mathnet format
            MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix a = new MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix(a1.Dimensions[0], a1.Dimensions[1]);
            for (int r = 0; r < a1.Dimensions[0]; r++)
                for (int c = 0; c < a1.Dimensions[1]; c++)
                    a[r, c] = new System.Numerics.Complex(
                        a1.InternalArray4Experts[c * a1.Dimensions[0] + r].real,
                        a1.InternalArray4Experts[c * a1.Dimensions[0] + r].imag);
            MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix b = new MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix(a2.Dimensions[0], a2.Dimensions[1]);
            for (int r = 0; r < a2.Dimensions[0]; r++)
                for (int c = 0; c < a2.Dimensions[1]; c++)
                    b[r, c] = new System.Numerics.Complex(
                        a2.InternalArray4Experts[c * a2.Dimensions[0] + r].real,
                        a2.InternalArray4Experts[c * a2.Dimensions[0] + r].imag);

            // Solve using QR decomposition
            int numQRs = b.ColumnCount;
            MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix x = new MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix(b.RowCount, b.ColumnCount);
            // Ax = b  =>  Rx = Q'b
            MathNet.Numerics.LinearAlgebra.Complex.Factorization.DenseQR qr = new MathNet.Numerics.LinearAlgebra.Complex.Factorization.DenseQR(a);
            MathNet.Numerics.LinearAlgebra.Generic.Matrix<System.Numerics.Complex> qt = qr.Q.ConjugateTranspose();
            int n = qr.R.ColumnCount;
            for (int i = 0; i < numQRs; i++)
            {
                // bi = b(:,i)
                MathNet.Numerics.LinearAlgebra.Generic.Matrix<System.Numerics.Complex> bi = b.SubMatrix(0, b.RowCount, i, 1);
                // Rx = c (where c = Q'b)
                MathNet.Numerics.LinearAlgebra.Generic.Matrix<System.Numerics.Complex> qb = qt * bi;
                // Back substitute
                for (int rj = 0; rj < x.RowCount; rj++)
                {
                    int j = x.RowCount - rj - 1;
                    for (int k = j; k < n; k++)
                        x[j, i] += qr.R[j, k] * x[k, i];
                    x[j, i] = (qb[j, 0] - x[j, i]) / qr.R[j, j];
                }
            }

            ILArray<complex> outArray = new ILArray<complex>(a2.Dimensions.ToIntArray());

            for (int r = 0; r < outArray.Dimensions[0]; r++)
                for (int c = 0; c < outArray.Dimensions[1]; c++)
                    outArray.InternalArray4Experts[c * outArray.Dimensions[0] + r] = new complex(
                        x[r, c].Real,
                        x[r, c].Imaginary);

            return outArray;
        }

        /// <summary>
        /// Performs left division of the two supplied complex arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Left-division of two arrays</returns>
        public static ILArray<fcomplex> ldivide(ILArray<fcomplex> a1, ILArray<fcomplex> a2)
        {
            // Convert matrices to mathnet format
            MathNet.Numerics.LinearAlgebra.Complex32.DenseMatrix a = new MathNet.Numerics.LinearAlgebra.Complex32.DenseMatrix(a1.Dimensions[0], a1.Dimensions[1]);
            for (int r = 0; r < a1.Dimensions[0]; r++)
                for (int c = 0; c < a1.Dimensions[1]; c++)
                    a[r, c] = new MathNet.Numerics.Complex32(
                        a1.InternalArray4Experts[c * a1.Dimensions[0] + r].real,
                        a1.InternalArray4Experts[c * a1.Dimensions[0] + r].imag);
            MathNet.Numerics.LinearAlgebra.Complex32.DenseMatrix b = new MathNet.Numerics.LinearAlgebra.Complex32.DenseMatrix(a2.Dimensions[0], a2.Dimensions[1]);
            for (int r = 0; r < a2.Dimensions[0]; r++)
                for (int c = 0; c < a2.Dimensions[1]; c++)
                    b[r, c] = new MathNet.Numerics.Complex32(
                        a2.InternalArray4Experts[c * a2.Dimensions[0] + r].real,
                        a2.InternalArray4Experts[c * a2.Dimensions[0] + r].imag);

            // Solve using QR decomposition
            int numQRs = b.ColumnCount;
            MathNet.Numerics.LinearAlgebra.Complex32.DenseMatrix x = new MathNet.Numerics.LinearAlgebra.Complex32.DenseMatrix(b.RowCount, b.ColumnCount);
            // Ax = b  =>  Rx = Q'b
            MathNet.Numerics.LinearAlgebra.Complex32.Factorization.DenseQR qr = new MathNet.Numerics.LinearAlgebra.Complex32.Factorization.DenseQR(a);
            MathNet.Numerics.LinearAlgebra.Generic.Matrix<MathNet.Numerics.Complex32> qt = qr.Q.ConjugateTranspose();
            int n = qr.R.ColumnCount;
            for (int i = 0; i < numQRs; i++)
            {
                // bi = b(:,i)
                MathNet.Numerics.LinearAlgebra.Generic.Matrix<MathNet.Numerics.Complex32> bi = b.SubMatrix(0, b.RowCount, i, 1);
                // Rx = c (where c = Q'b)
                MathNet.Numerics.LinearAlgebra.Generic.Matrix<MathNet.Numerics.Complex32> qb = qt * bi;
                // Back substitute
                for (int rj = 0; rj < x.RowCount; rj++)
                {
                    int j = x.RowCount - rj - 1;
                    for (int k = j; k < n; k++)
                        x[j, i] += qr.R[j, k] * x[k, i];
                    x[j, i] = (qb[j, 0] - x[j, i]) / qr.R[j, j];
                }
            }

            ILArray<fcomplex> outArray = new ILArray<fcomplex>(a2.Dimensions.ToIntArray());

            for (int r = 0; r < outArray.Dimensions[0]; r++)
                for (int c = 0; c < outArray.Dimensions[1]; c++)
                    outArray.InternalArray4Experts[c * outArray.Dimensions[0] + r] = new fcomplex(
                        x[r, c].Real,
                        x[r, c].Imaginary);

            return outArray;
        }

        /// <summary>
        /// Performs left division of the two supplied complex arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Left-division of two arrays</returns>
        public static ILArray<double> ldivide(ILArray<double> a1, ILArray<double> a2)
        {
            // Convert matrices to mathnet format
            MathNet.Numerics.LinearAlgebra.Double.DenseMatrix a = new MathNet.Numerics.LinearAlgebra.Double.DenseMatrix(a1.Dimensions[0], a1.Dimensions[1]);
            for (int r = 0; r < a1.Dimensions[0]; r++)
                for (int c = 0; c < a1.Dimensions[1]; c++)
                    a[r, c] = a1.InternalArray4Experts[c * a1.Dimensions[0] + r];
            MathNet.Numerics.LinearAlgebra.Double.DenseMatrix b = new MathNet.Numerics.LinearAlgebra.Double.DenseMatrix(a2.Dimensions[0], a2.Dimensions[1]);
            for (int r = 0; r < a2.Dimensions[0]; r++)
                for (int c = 0; c < a2.Dimensions[1]; c++)
                    b[r, c] = a2.InternalArray4Experts[c * a2.Dimensions[0] + r];

            // Solve using QR decomposition
            int numQRs = b.ColumnCount;
            MathNet.Numerics.LinearAlgebra.Double.DenseMatrix x = new MathNet.Numerics.LinearAlgebra.Double.DenseMatrix(b.RowCount, b.ColumnCount);
            // Ax = b  =>  Rx = Q'b
            MathNet.Numerics.LinearAlgebra.Double.Factorization.DenseQR qr = new MathNet.Numerics.LinearAlgebra.Double.Factorization.DenseQR(a);
            MathNet.Numerics.LinearAlgebra.Generic.Matrix<double> qt = qr.Q.ConjugateTranspose();
            int n = qr.R.ColumnCount;
            for (int i = 0; i < numQRs; i++)
            {
                // bi = b(:,i)
                MathNet.Numerics.LinearAlgebra.Generic.Matrix<double> bi = b.SubMatrix(0, b.RowCount, i, 1);
                // Rx = c (where c = Q'b)
                MathNet.Numerics.LinearAlgebra.Generic.Matrix<double> qb = qt * bi;
                // Back substitute
                for (int rj = 0; rj < x.RowCount; rj++)
                {
                    int j = x.RowCount - rj - 1;
                    for (int k = j; k < n; k++)
                        x[j, i] += qr.R[j, k] * x[k, i];
                    x[j, i] = (qb[j, 0] - x[j, i]) / qr.R[j, j];
                }
            }

            ILArray<double> outArray = new ILArray<double>(a2.Dimensions.ToIntArray());

            for (int r = 0; r < outArray.Dimensions[0]; r++)
                for (int c = 0; c < outArray.Dimensions[1]; c++)
                    outArray.InternalArray4Experts[c * outArray.Dimensions[0] + r] = x[r, c];
            ILBaseArray<double> zz = a1.ShiftDimensions(1);
            return outArray;
        }

        /// <summary>
        /// Performs left division of the two supplied complex arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Left-division of two arrays</returns>
        public static ILArray<float> ldivide(ILArray<float> a1, ILArray<float> a2)
        {
            // Convert matrices to mathnet format
            MathNet.Numerics.LinearAlgebra.Single.DenseMatrix a = new MathNet.Numerics.LinearAlgebra.Single.DenseMatrix(a1.Dimensions[0], a1.Dimensions[1]);
            for (int r = 0; r < a1.Dimensions[0]; r++)
                for (int c = 0; c < a1.Dimensions[1]; c++)
                    a[r, c] = a1.InternalArray4Experts[c * a1.Dimensions[0] + r];
            MathNet.Numerics.LinearAlgebra.Single.DenseMatrix b = new MathNet.Numerics.LinearAlgebra.Single.DenseMatrix(a2.Dimensions[0], a2.Dimensions[1]);
            for (int r = 0; r < a2.Dimensions[0]; r++)
                for (int c = 0; c < a2.Dimensions[1]; c++)
                    b[r, c] = a2.InternalArray4Experts[c * a2.Dimensions[0] + r];

            // Solve using QR decomposition
            int numQRs = b.ColumnCount;
            MathNet.Numerics.LinearAlgebra.Single.DenseMatrix x = new MathNet.Numerics.LinearAlgebra.Single.DenseMatrix(b.RowCount, b.ColumnCount);
            // Ax = b  =>  Rx = Q'b
            MathNet.Numerics.LinearAlgebra.Single.Factorization.DenseQR qr = new MathNet.Numerics.LinearAlgebra.Single.Factorization.DenseQR(a);
            MathNet.Numerics.LinearAlgebra.Generic.Matrix<float> qt = qr.Q.ConjugateTranspose();
            int n = qr.R.ColumnCount;
            for (int i = 0; i < numQRs; i++)
            {
                // bi = b(:,i)
                MathNet.Numerics.LinearAlgebra.Generic.Matrix<float> bi = b.SubMatrix(0, b.RowCount, i, 1);
                // Rx = c (where c = Q'b)
                MathNet.Numerics.LinearAlgebra.Generic.Matrix<float> qb = qt * bi;
                // Back substitute
                for (int rj = 0; rj < x.RowCount; rj++)
                {
                    int j = x.RowCount - rj - 1;
                    for (int k = j; k < n; k++)
                        x[j, i] += qr.R[j, k] * x[k, i];
                    x[j, i] = (qb[j, 0] - x[j, i]) / qr.R[j, j];
                }
            }

            ILArray<float> outArray = new ILArray<float>(a2.Dimensions.ToIntArray());

            for (int r = 0; r < outArray.Dimensions[0]; r++)
                for (int c = 0; c < outArray.Dimensions[1]; c++)
                    outArray.InternalArray4Experts[c * outArray.Dimensions[0] + r] = x[r, c];

            return outArray;
        }

        /// <summary>
        /// Performs left division of the two supplied complex arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Left-division of two arrays</returns>
        public static ILArray<double> ldivide(ILArray<long> a1, ILArray<long> a2)
        {
            return ldivide(ILMath.todouble(a1), ILMath.todouble(a2));
        }

        /// <summary>
        /// Performs left division of the two supplied complex arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Left-division of two arrays</returns>
        public static ILArray<double> ldivide(ILArray<ulong> a1, ILArray<ulong> a2)
        {
            return ldivide(ILMath.todouble(a1), ILMath.todouble(a2));
        }

        /// <summary>
        /// Performs left division of the two supplied complex arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Left-division of two arrays</returns>
        public static ILArray<double> ldivide(ILArray<int> a1, ILArray<int> a2)
        {
            return ldivide(ILMath.todouble(a1), ILMath.todouble(a2));
        }

        /// <summary>
        /// Performs left division of the two supplied complex arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Left-division of two arrays</returns>
        public static ILArray<double> ldivide(ILArray<uint> a1, ILArray<uint> a2)
        {
            return ldivide(ILMath.todouble(a1), ILMath.todouble(a2));
        }

        /// <summary>
        /// Performs left division of the two supplied complex arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Left-division of two arrays</returns>
        public static ILArray<double> ldivide(ILArray<short> a1, ILArray<short> a2)
        {
            return ldivide(ILMath.todouble(a1), ILMath.todouble(a2));
        }

        /// <summary>
        /// Performs left division of the two supplied complex arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Left-division of two arrays</returns>
        public static ILArray<double> ldivide(ILArray<ushort> a1, ILArray<ushort> a2)
        {
            return ldivide(ILMath.todouble(a1), ILMath.todouble(a2));
        }

        /// <summary>
        /// Performs left division of the two supplied complex arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Left-division of two arrays</returns>
        public static ILArray<double> ldivide(ILArray<char> a1, ILArray<char> a2)
        {
            return ldivide(ILMath.todouble(a1), ILMath.todouble(a2));
        }

        /// <summary>
        /// Performs left division of the two supplied complex arrays.
        /// </summary>
        /// <param name="a1">First array</param>
        /// <param name="a2">Second array</param>
        /// <returns>Left-division of two arrays</returns>
        public static ILArray<double> ldivide(ILArray<byte> a1, ILArray<byte> a2)
        {
            return ldivide(ILMath.todouble(a1), ILMath.todouble(a2));
        }

        public static void insertMatrix<T>(ILArray<T> target, IEnumerable<IronPython.Runtime.Slice> ranges, ILArray<T> source)
        {
            int tgtCols = target.Dimensions[0];
            int srcCols = source.Dimensions[0];

            int[] starts = new int[2];
            int[] stops = new int[2];
            int i = 0;
            foreach (IronPython.Runtime.Slice range in ranges)
            {
                starts[i] = (int)range.start;
                stops[i] = (int)range.stop;
                i++;
            }

            for (int j = starts[1]; j < stops[1]; j++)
                for (i = starts[0]; i < stops[0]; i++)
                    target.InternalArray4Experts[j * tgtCols + i] = source.InternalArray4Experts[(j - starts[1]) * srcCols + (i - starts[0])];
        }

        public static void insertMatrix(ILArray<double> target, IEnumerable<int> position, double source)
        {
            int[] pos = new int[2];
            int i = 0;
            foreach (int p in position) pos[i++] = p;

            target.InternalArray4Experts[pos[1] * target.Dimensions[0] + pos[0]] = source;
        }

        public static void insertMatrix(ILArray<complex> target, IEnumerable<int> position, complex source)
        {
            int[] pos = new int[2];
            int i = 0;
            foreach (int p in position) pos[i++] = p;

            target.InternalArray4Experts[pos[1] * target.Dimensions[0] + pos[0]] = source;
        }

        public static void insertMatrix(ILArray<float> target, IEnumerable<int> position, float source)
        {
            int[] pos = new int[2];
            int i = 0;
            foreach (int p in position) pos[i++] = p;

            target.InternalArray4Experts[pos[1] * target.Dimensions[0] + pos[0]] = source;
        }

        public static void insertMatrix(ILArray<fcomplex> target, IEnumerable<int> position, fcomplex source)
        {
            int[] pos = new int[2];
            int i = 0;
            foreach (int p in position) pos[i++] = p;

            target.InternalArray4Experts[pos[1] * target.Dimensions[0] + pos[0]] = source;
        }

        public static void insertMatrix(ILArray<int> target, IEnumerable<int> position, int source)
        {
            int[] pos = new int[2];
            int i = 0;
            foreach (int p in position) pos[i++] = p;

            target.InternalArray4Experts[pos[1] * target.Dimensions[0] + pos[0]] = source;
        }

        public static void insertMatrix(ILArray<uint> target, IEnumerable<int> position, uint source)
        {
            int[] pos = new int[2];
            int i = 0;
            foreach (int p in position) pos[i++] = p;

            target.InternalArray4Experts[pos[1] * target.Dimensions[0] + pos[0]] = source;
        }

        public static void insertMatrix(ILArray<byte> target, IEnumerable<int> position, byte source)
        {
            int[] pos = new int[2];
            int i = 0;
            foreach (int p in position) pos[i++] = p;

            target.InternalArray4Experts[pos[0] * target.Dimensions[1] + pos[1]] = source;
        }

        public static void insertMatrix(ILArray<char> target, IEnumerable<int> position, char source)
        {
            int[] pos = new int[2];
            int i = 0;
            foreach (int p in position) pos[i++] = p;

            target.InternalArray4Experts[pos[1] * target.Dimensions[0] + pos[0]] = source;
        }

        public static void insertMatrix(ILArray<bool> target, IEnumerable<int> position, bool source)
        {
            int[] pos = new int[2];
            int i = 0;
            foreach (int p in position) pos[i++] = p;

            target.InternalArray4Experts[pos[1] * target.Dimensions[0] + pos[0]] = source;
        }

        public static ILArray<double> expandArray(ILArray<double> array, IEnumerable<int> newDims)
        {
            // Convert newDims to int[]
            int numDims = 0;
            foreach (int dim in newDims) numDims++;
            int[] newDimsArray = new int[numDims];
            int d = 0;
            foreach (int dim in newDims) newDimsArray[d++] = dim;

            // Create new array
            ILArray<double> newArray = new ILArray<double>(newDimsArray);

            // Copy from old array into new array
            for (int i = 0; i < array.InternalArray4Experts.Length; i++)
            {
                int[] inSub = ILMath.ind2sub(array, i);
                int outInd = ILMath.sub2ind((ILArray<double>)newArray, inSub);
                newArray.InternalArray4Experts[outInd] = array.InternalArray4Experts[i];
            }

            // Return expanded array
            return newArray;
        }

        /// <summary>
        /// Returns a specified slice of an array.
        /// </summary>
        /// <param name="array">Array to slice</param>
        /// <param name="indices">Indices describing slice to return</param>
        /// <returns>Sliced array</returns>
        public static ILArray<T> getArraySlice<T>(ILArray<T> array, IronPython.Runtime.List indices)
        {
            int dest = 0;
            ILArray<T> outArray;

            if (array.Dimensions.NonSingletonDimensions > 1 && indices.Count == 1)
            {
                // Flat indexing into multi-dimensional array
                int start, step, stop, size;
                if (indices[0] is IronPython.Runtime.List)
                {
                    // Range defined by a slice
                    IronPython.Runtime.List slice = indices[0] as IronPython.Runtime.List;
                    start = (int)double.Parse(slice[0].ToString()) - 1;
                    step = (int)double.Parse(slice[1].ToString());
                    double stopVal;
                    if (double.TryParse(slice[2].ToString(), out stopVal))
                    {
                        // Bounded slice
                        stop = (int)stopVal;
                    }
                    else
                    {
                        // Unbounded slice
                        stop = array.Dimensions.NumberOfElements - 1;
                    }
                    size = (stop - start) / step + 1;
                }
                /*else if (indices[0] is ILArray<double>)
                {
                    ILArray<double> indexArr = indices[0] as ILArray<double>;
                    size = indexArr.Length;
                    start = (int)indexArr[0] - 1;
                    if (size > 1)
                        stop = (int)indexArr[size - 1];
                    else
                        stop = start;
                    step = (stop - start + 1) / size;
                }*/
                else
                {
                    // Single value
                    start = (int)double.Parse(indices[0].ToString());
                    step = 1;
                    stop = start + 1;
                    size = 1;
                }

                outArray = new ILArray<T>(size);
                if (step > 0 && stop >= start)
                    for (int n = 0, m = start; m <= stop; m += step, n++)
                        outArray.InternalArray4Experts[n] = array.InternalArray4Experts[m];
                else if (step < 0 && stop <= start)
                    for (int n = 0, m = stop; m >= start; m += step, n++)
                        outArray.InternalArray4Experts[n] = array.InternalArray4Experts[m];
            }
            else
            {
                int N = array.Dimensions.NumberOfDimensions;
                int[] starts = new int[N];
                int[] steps = new int[N];
                int[] stops = new int[N];
                int[] sizes = new int[N];

                // Regular indexing
                int curIdx = 0;
                for (int i = 0; i < N; i++)
                {
                    if (array.Dimensions[i] == 1)
                    {
                        // Singleton dimension
                        starts[i] = 0;
                        steps[i] = 1;
                        stops[i] = 0;
                        sizes[i] = 1;
                    }
                    else
                    {
                        if (indices[curIdx] is IronPython.Runtime.List)
                        {
                            // Range defined by a slice
                            IronPython.Runtime.List slice = indices[curIdx] as IronPython.Runtime.List;
                            starts[i] = (int)double.Parse(slice[0].ToString()) - 1;
                            steps[i] = (int)double.Parse(slice[1].ToString());
                            double stopVal;
                            if (double.TryParse(slice[2].ToString(), out stopVal))
                            {
                                // Bounded slice
                                stops[i] = (int)stopVal - 1;
                            }
                            else
                            {
                                // Unbounded slice
                                stops[i] = array.Dimensions[i] - 1;
                            }
                            sizes[i] = (stops[i] - starts[i]) / steps[i] + 1;
                        }
                        /*else if (indices[curIdx] is ILArray<double>)
                        {
                            ILArray<double> indexArr = indices[curIdx] as ILArray<double>;
                            sizes[i] = indexArr.Length;
                            starts[i] = (int)indexArr[0] - 1;
                            if (sizes[i] > 1)
                                stops[i] = (int)indexArr[sizes[i] - 1];
                            else
                                stops[i] = starts[i];
                            steps[i] = (stops[i] - starts[i] + 1) / sizes[i];
                        }*/
                        else
                        {
                            // Single value
                            stops[i] = (int)double.Parse(indices[curIdx].ToString()) - 1;
                            steps[i] = 1;
                            starts[i] = stops[i];
                            sizes[i] = 1;
                        }

                        curIdx++;
                    }
                }

                // Create sub-array
                outArray = new ILArray<T>(sizes);

                // Set up indices
                int[] idx = new int[N];
                for (int d = 0; d < N; d++)
                    idx[d] = starts[d];

                int[] dims = array.Dimensions.ToIntArray();

                // Copy data into sub-array               
                while (dest < outArray.InternalArray4Experts.Length)
                {
                    // Calculate flat source index
                    int n = idx.Length;
                    int src = idx[n - 1];

                    int nd = n - 2;

                    while (nd > -1)
                    {
                        src *= dims[nd];
                        src += idx[nd];
                        nd--;
                    }

                    outArray.InternalArray4Experts[dest] = array.InternalArray4Experts[src];

                    n = idx.Length;

                    int d = 0;

                    while ( d < n)
                    {
                        idx[d] += steps[d];
                        if (idx[d] <= stops[d]) break;

                        idx[d] = starts[d];
                        d++;
                    }

                    dest++;
                }
            }

            return outArray;
        }

        /// <summary>
        /// Set the contents of an array slice.
        /// </summary>
        /// <param name="destn">Destination array to receive source array</param>
        /// <param name="indices">Indices describing slice</param>
        /// <param name="src">Source array to set into the destination array</param>
        public static void setArraySlice<T>(ILArray<T> destn, IronPython.Runtime.List indices, ILArray<T> src)
        {
            if (destn.Dimensions.NonSingletonDimensions > 1 && indices.Count == 1)
            {
                // Flat indexing into multi-dimensional array
                int start, step, stop, size;
                if (indices[0] is IronPython.Runtime.List)
                {
                    // Range defined by a slice
                    IronPython.Runtime.List slice = indices[0] as IronPython.Runtime.List;
                    start = (int)double.Parse(slice[0].ToString()) - 1;
                    step = (int)double.Parse(slice[1].ToString());
                    double stopVal;
                    if (double.TryParse(slice[2].ToString(), out stopVal))
                    {
                        // Bounded slice
                        stop = (int)stopVal - 1;
                    }
                    else
                    {
                        // Unbounded slice
                        stop = destn.Dimensions.NumberOfElements - 1;
                    }
                    size = (stop - start) / step + 1;
                }
                else if (indices[0] is ILArray<double>)
                {
                    ILArray<double> indexArr = indices[0] as ILArray<double>;
                    size = indexArr.Length;
                    start = (int)indexArr[0] - 1;
                    if (size > 1)
                        stop = (int)indexArr[size - 1];
                    else
                        stop = start;
                    step = (stop - start + 1) / size;
                }
                else
                {
                    // Single value
                    start = (int)double.Parse(indices[0].ToString());
                    step = 1;
                    stop = start + 1;
                    size = 1;
                }

                if (step > 0 && stop >= start)
                    for (int n = 0, m = start; m <= stop; m += step, n++)
                        destn.InternalArray4Experts[m] = src.InternalArray4Experts[n];
                else if (step < 0 && stop <= start)
                    for (int n = 0, m = stop; m >= start; m += step, n++)
                        destn.InternalArray4Experts[m] = src.InternalArray4Experts[n];
            }
            else
            {
                int N = destn.Dimensions.NumberOfDimensions;
                int[] starts = new int[N];
                int[] steps = new int[N];
                int[] stops = new int[N];
                int[] sizes = new int[N];

                // Regular indexing
                int curIdx = 0;
                for (int i = 0; i < N; i++)
                {
                    if (destn.Dimensions[i] == 1)
                    {
                        // Singleton dimension
                        starts[i] = 0;
                        steps[i] = 1;
                        stops[i] = 0;
                        sizes[i] = 1;
                    }
                    else
                    {
                        if (indices[curIdx] is IronPython.Runtime.List)
                        {
                            // Range defined by a slice
                            IronPython.Runtime.List slice = indices[curIdx] as IronPython.Runtime.List;
                            starts[i] = (int)double.Parse(slice[0].ToString()) - 1;
                            steps[i] = (int)double.Parse(slice[1].ToString());
                            double stopVal;
                            if (double.TryParse(slice[2].ToString(), out stopVal))
                            {
                                // Bounded slice
                                stops[i] = (int)stopVal;
                            }
                            else
                            {
                                // Unbounded slice
                                stops[i] = destn.Dimensions[i] - 1;
                            }
                            sizes[i] = (stops[i] - starts[i]) / steps[i] + 1;
                        }
                        /*else if (indices[curIdx] is ILArray<double>)
                        {
                            ILArray<double> indexArr = indices[curIdx] as ILArray<double>;
                            sizes[i] = indexArr.Length;
                            starts[i] = (int)indexArr[0] - 1;
                            if (sizes[i] > 1)
                                stops[i] = (int)indexArr[sizes[i] - 1];
                            else
                                stops[i] = starts[i];
                            steps[i] = (stops[i] - starts[i] + 1) / sizes[i];
                        }*/
                        else
                        {
                            // Single value
                            stops[i] = (int)double.Parse(indices[curIdx].ToString()) - 1;
                            steps[i] = 1;
                            starts[i] = stops[i];
                            sizes[i] = 1;
                        }

                        curIdx++;
                    }
                }

                // Set up indices
                int[] idx = new int[N];
                for (int d = 0; d < N; d++)
                    idx[d] = starts[d];

                int[] dims = destn.Dimensions.ToIntArray();

                // Copy data into sub-array
                for (int srcIdx = 0; srcIdx < src.InternalArray4Experts.Length; srcIdx++)
                {
                    // Calculate flat source index
                    int destnIdx = Utils.indicesToFlat(dims, idx);
                    destn.InternalArray4Experts[destnIdx] = src.InternalArray4Experts[srcIdx];
                    Utils.indicesInc(starts, steps, stops, ref idx);
                }
            }
        }

        /// <summary>
        /// Transpose matrix.
        /// </summary>
        /// <typeparam name="T">Element type of matrix</typeparam>
        /// <param name="a1">Matrix to transpose</param>
        /// <param name="conjugate">Whether to conjugate complex numbers</param>
        /// <returns>Transposed matrix</returns>
        public static ILBaseArray<T> transpose<T>(ILBaseArray<T> a1, bool conjugate)
        {
            ILBaseArray<T> outArray = new ILArray<T>(a1.Dimensions[1], a1.Dimensions[0]);

            for (int r = 0; r < a1.Dimensions[0]; r++)
                for (int c = 0; c < a1.Dimensions[1]; c++)
                    outArray.InternalArray4Experts[r * a1.Dimensions[1] + c] = a1.InternalArray4Experts[c * a1.Dimensions[0] + r];

            if (conjugate)
            {
                if (outArray is ILArray<complex>)
                {
                    ILBaseArray<complex> o = outArray as ILBaseArray<complex>;
                    for (int r = 0; r < a1.Dimensions[0]; r++)
                        for (int c = 0; c < a1.Dimensions[1]; c++)
                            o.InternalArray4Experts[r * a1.Dimensions[1] + c].imag = -o.InternalArray4Experts[r * a1.Dimensions[1] + c].imag;
                }
                else if (outArray is ILArray<fcomplex>)
                {
                    ILBaseArray<fcomplex> o = outArray as ILBaseArray<fcomplex>;
                    for (int r = 0; r < a1.Dimensions[0]; r++)
                        for (int c = 0; c < a1.Dimensions[1]; c++)
                            o.InternalArray4Experts[r * a1.Dimensions[1] + c].imag = -o.InternalArray4Experts[r * a1.Dimensions[1] + c].imag;
                }
            }

            return outArray;
        }
    }
}
